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Chapter 1 
Introduction 


This is the final report on NASA Grant No. NAG 1-1479 entitled, “Compu- 
tational Aeroacoustics and Numerical Simulation of Supersonic Jets.” The 
Principal Investigators are Drs. Philip J. Morris and Lyle N. Long of the 
Department of Aerospace Engineering at the Pennsylvania State University. 

The research project has been a computational study of computational 
aeroacoustics algorithms and numerical simulations of the flow and noise of 
supersonic jets. Dining this study a new method for the implementation 
of solid wall boundary conditions for complex geometries in three dimen- 
sions has been developed. In addition, a detailed study of the simulation 
of the flow in and noise from supersonic circular and rectangular jets has 
been conducted. Extensive comparisons have been made with experimental 
measurements. A summary of the results of the research program are at- 
tached as the main body of this report in the form of two publications. In 
addition, the report lists the names of the students who were supported by 
this grant, their degrees, and the titles of their dissertations. In addition, 
a list of presentations and publications made by the Principal Investigators 
and the research students is also included. 
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Chapter 2 


Publications and Presentations 


Morris, P. J., Long, L. N., Chung, C. and Chyczewski, T., “Computa- 
tional aeroacoustics algorithms: nonuniform grids,” AIAA Paper 94-2295, 
25th AIAA Fluid Dynamics Conference, June 1994. 

Chung, C. and Morris, P. J., “Wave propagation and scattering in com- 
putational aeroacoustics,” ICASE/LaRC Workshop in Benchmark Problems 
in Computational Aeroacoustics, Hampton, VA, October 1994. 

Morris, P. J., Chung, C. and Pautet, L. R., “Acoustic scattering: a nu- 
merical simulation,” 47th APS/Division of Fluid Dynamics Meeting, Atlanta, 
GA, November 1994. 

Chung, C. and Morris, P. J., “A new boundary treatment for two- and 
three-dimensional acoustic scattering problems,” AIAA/CEAS Aeroacoustics 
Conference, Munich, Germany, June 1995. 

Chung, C. and Morris, P. J., “Acoustic scattering from two- and three- 
dimensional bodies,” submitted for publication to J. Computational Acous- 
tics, August 1996. (copy attached below) 

Chyczewski, T. and Long, L. N., “Numerical prediction of the noise pro- 
duced by a perfectly-expanded rectangular jet,” 2nd AIAA/CEAS Aeroa- 
coustics Conference, State College, PA, May 1996. (copy attached below) 
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Chapter 3 

Graduate Student Researchers 


Tom Chyczewski, “A time-dependent, three-dimensional numerical study 
of supersonic rectangular jet flow and noise using the full Navier-Stokes equa- 
tions,” Ph. D. thesis, Aerospace Engineering, 1996. Present Occupation: 
Research Engineer at Allison Gas Turbines, Indianapolis, IN. 

Cathy Chung, “Wave propagation and scattering in computational aeroa- 
coustics,” Ph.D. thesis, Aerospace Engineering, 1995. Present Occupation: 
Homemaker and mother. 
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Appendix: Summary Papers 
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In this paper we consider the scattering of sound by two- and three-dimensional bodies with 
arbitrary geometries. Particular emphasis is placed on the methodology for the implementation of 
solid wall boundary conditions for high-order, high-bandwidth numerical schemes. The Impedance 
Mismatch Method (IMM) is introduced to treat solid wall boundaries. In this method the solid 
wall is simulated using a wall region in which the characteristic impedance is set to a different 
value from that in the fluid region. This method has many advantages over traditional solid wall 
boundary treatments, including simplicity of coding, speed of computation and the ability to treat 
curved boundaries. This method has been used to solve a number of acoustic scattering problems 
to demonstrate its effectiveness. These problems include acoustic reflections from an infinite plate, 
acoustic scattering from a two-dimensional finite plate and a cylinder, and acoustic scattering by 
a sphere and a cylindrical shell. 


1. Introduction 

Acoustic scattering by solid bodies is of importance in engineering noise prediction and 
control. In order to simulate acoustic scattering accurately, we need high-order numerical 
schemes as well as methods to implement solid wall boundary conditions for such high-order 
schemes. In recent years, a number of high-order schemes have been developed. One of them 
is the Dispersion- Relation-Preserving (DRP) scheme developed by Tam and Webb 7 . This 
scheme is used in this paper. 

The most important consideration in computations of acoustic scattering by solid bodies 
is the implementation of solid wall boundary conditions. For inviscid flow the solid wall 
boundary condition requires that the velocity normal to the wall is zero. For low-order finite 
difference schemes or finite volume schemes, the imposition of solid wall boundary conditions 
can usually be carried out in a straightforward manner (Khan et al 3 , Huh et al 2 );Though, 
it should be noted that in these references, the problem was reformulated to solve for the 
scattered field only. For high-order finite difference schemes, treatment of the wall condition 
is complicated and it has had little investigation. Recently, Tam and Dong 6 proposed a way 
to implement solid wall boundary conditions for high-order finite difference DRP schemes. 
In their method, ghost points are needed and the seven-point central difference spatial stencil 
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must be changed to a one-sided stencil when computations are performed at grid points 
within three rows of the solid surface. For an object with a simple geometry, such as a flat 
plate, this solid wall boundary condition can be used relatively easily. But for a complicated 
surface where the boundary is curved, the implementation of this wall boundary condition 
involves considerable work in coding and long computer times. This is especially true if 
the calculations are to be performed on computers with parallel architecture. Kurbatskii 
and Tam 4 have extended the plane wall boundary method of Tam and Dong 6 to two- 
dimensional curved walls. 

In this paper we introduce a very efficient method to implement solid wall boundary 
conditions. This is the Impedance Mismatch Method (IMM). This method can be applied 
easily to high-order finite difference schemes. In this method the solid wall is simulated 
using a wall region in which the characteristic impedance is set to a different value from 
that in the fluid region. When acoustic waves encounter the interface between the two 
different regions, they will be reflected in-phase; the interface acts just like a solid surface. 
Actually, in this method, no wall boundary conditions need to be implemented, all that 
is needed is to define a body region and to set a different characteristic impedance in this 
region. Since this method does not involve any changes in the stencil, it can be used to 
represent the geometry of any object without difficulty. Also it makes the computation 
much faster and coding much simpler compared to traditional ways of dealing with solid 
wall boundaries. 

In the next section of this paper, the one-dimensional linearized Euler equations are used 
to introduce the Impedance Mismatch Method. Then several numerical simulations are 
performed for acoustic scattering by two-dimensional bodies, including acoustic reflection 
from an infinite plate, and acoustic scattering from a finite plate and from a cylinder. 
Following that, acoustic scattering by three-dimensional bodies is computed. Examples 
include acoustic scattering by a sphere and a cylindrical shell. The numerical results are 
compared with either analytical solutions or solutions obtained using traditional solid wall 
boundary conditions. Finally, the advantages and some disadvantages of the IMM are given. 

2. The Impedance Mismatch Method (IMM) 

In this section a one-dimensional example is used to fix the idea of the IMM. First, the 
elementary problem of plane wave reflections is presented to introduce the characteristic 
impedance. However, it is found that when a direct simulation of this problem is extended 
to two dimensions a numerical instability occurs. In order to overcome this difficulty, an 
auxiliary problem is proposed. Analysis shows that the auxiliary problem gives the same 
solution as the physical problem in the region of interest external to the body. 

2.1. The Physical Problem 

From classical acoustics theory, it is known that when a plane wave in a fluid medium 
impinges normally to the boundary of a contiguous second medium, a reflected wave is 
generated in the first medium and a transmitted wave moves into the second medium. The 
ratio of the pressure amplitude of the reflected wave to that of incident wave depends on 
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the characteristic impedances, pa, 
by 


of the two media, the pressure amplitude ratio is given 


\Pr/Pi\ 


P2&2 ~ Pl°l 
P2 a 2 + Pl a l 


( 2 . 1 ) 


where pi, p r are the incident and reflected pressures respectively; pi and P 2 are the mean 
densities in the first and second media respectively, and oi and 02 are respective speeds 
of sound. This relation is obtained using the conditions that the pressure and the particle 
velocity are continuous at the boundary between the two media. This means that the 
normal derivatives of the pressure and the particle velocity are discontinuous. When the 
second medium has a much higher characteristic impedance, most of the wave energy is 
reflected. As the ratio of the characteristic impedance of the second medium to the first 
approaches infinity, all the incident waves are reflected. The second medium acts like a solid 
object. Thus, setting a higher impedance in a certain region can be used to simulate the 
effect of a solid object in this region. This is the basic idea behind the IMM. In the present 
formulation the speed of sound inside and outside the body is kept the same. This means 
that the wave speed is constant throughout the domain and permits the CFL number to be 
kept at almost the same value as when no object is present. Thus, it is the mean density 
in the wall region that is modified to provide the impedance mismatch. 

This problem is governed by the linearized Euler equations. The one-dimensional form 
of these equations without mean flow is 




du 1 d\ 
dt po d : 

dp 9 d' 

m +poa °d: 


( 2 . 2 ) 

(2.3) 

(2.4) 


In above equations, all quantities are nondimensionalized by characteristic scales: the mesh 
size of the numerical simulation (Ax = Ay) for a length scale; the ambient speed of sound 
(a amb ) for a velocity scale; the ambient density ( p am b ) for a density scale; pambO- 2 amb for a 
pressure scale, po is the nondimensionalized mean density and do is the nondimensionalized 
local mean sound speed, do is equal to unity at all points for a constant speed of sound and 
po is equal to unity exterior to the body. 

An initial value problem is now solved to test the IMM. The computational domain is 
-100 < x < 100. The wall region is 50 < x < 100, inside which the mean density, pi - 30. 
An acoustic pulse is generated in the center of the domain at t = 0. The initial conditions 


are 


t ln2 2x n 

p = p = exp(-— x ), u = 0 


(2.5) 


Figure 1 shows the computed pressure distribution at t = 75, At = 0.05. The analytical 
solution is also plotted. The first pulse from the left-travelling disturbance from the initial 
pulse. The second pulse is the reflected wave from the interface at x = 50, also traveling 
to the left. The third pulse is the transmitted wave, traveling to the right inside the wall 
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region. It can be seen that the pressure wave has been reflected after it encounters the 
interface at x = 50. The agreement between the computed and analytical solutions is good, 
but some small errors can still be seen in the plot. These errors are mainly due to the finite 
choice of the mean density in the wall region p 2 and ambiguities about the exact position 
of the wall boundary. These issues are discussed in later sections. 

The above example demonstrates that the IMM is applicable to this simple one- dimen- 
sional problem; however, when this method is used in multi-dimensional cases, instability 
occurs when the same time step is used. It is perhaps surprising that the one-dimensional 
problem gives any form of accurate solution. In Eqs. (2.3) and (2.4), the time derivatives 
are continuous functions. The coefficients and the spatial derivatives are only piecewise 
continuous and it is only in combination that they produce continuous functions. Finite- 
difference approximations depend on the smoothness of the function for accuracy. This is 
particularly true for high-order finite-difference approximations. Thus the finite-difference 
approximation to the spatial derivative will be inaccurate in the present case in the vicinity 
of the interface. So the large discontinuity in the coefficient (po/p a mb) will not t> e balanced 
by an accurate discontinuous behavior in the approximation to the spatial derivative. In 
order to overcome this problem, an auxiliary problem is introduced in which the discon- 
tinuous coefficients are combined with the primitive variables under the spatial derivative 
operator. 
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2.2. The Auxiliary Problem 

The 1-D linearized Euler equations without mean flow are re-written as 


dt dx 
du dp 

aT + fc-° 

dp du 

m + d^~° 


(2.6) 

(2.7) 

( 2 . 8 ) 


where p = p/po, u = upo, p = p/po- Note that ao equals unity. This set of equations is 
the same as that in the physical problem in the fluid region and wall region, but not at the 
interface, since po is only piecewise uniform and has a jump at interface. For this set of 
equations, the coefficients are continuous. The condition is imposed that variables p, u and 
p are also continuous at the interface of the two media. This is not the physical problem 
or interface condition. For a normal incident plane wave, the pressure amplitude of the 
reflected waves is found to be 


I Pr/Pi 


1/P2-1 
1/^2 + 1 


(2.9) 


If the wall region still mimics a solid wall in this case, then the mean density in the 
wall region must be set to a lower value instead of a higher value than that in the fluid 
region. This assumes that the mean speed of sound is kept the same in both the fluid and 
the wall regions. Since p 0 is always unity in the fluid region, the physical solution then can 
be obtained in this region from the solution of this auxiliary problem. 
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Fig. 2. Pressure at t = 75 for the 
one-dimensional auxiliary problem. 


Fig. 3. Pressure at t = 75 for the 
one-dimensional auxiliary problem for different 
density ratios. 


The same initial value problem has been solved by setting P 2 — 1/30 in the wall region. 
However, this time, the auxiliary problem is implemented. Figure 2 shows the computed 
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pressure distribution at t = 75. The analytical solution is also plotted. It can be seen 
that the physical solution is obtained in the fluid region, and the agreement between the 
computed and analytical solutions is good. 

Figure 3 also shows the pressure distributions for different mean density ratios in the wall 
region. It can be seen that the accuracy of computations depends on the density ratio; the 
smaller the ratio, the more closely the numerical simulation follows the rigid wall solution. 
But, it has been found numerically that this density ratio can not be set infinitely small 
due to the stability considerations. This is one of the sources of inaccuracies in the IMM. 


3. Governing Equations and Algorithms 

In this paper the governing equations to be solved are the linearized Euler equations. 
For three-dimensional flows with a mean flow only in x direction, they have the form 


where 


dU dE_ &F dG 
dt dx dy dz 
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0 
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0 
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P/P o 

. Pov . 


PqU) 


(3.10) 


M is the mean flow Mach number and H is a source term. All the quantities are non- 
dimensionalized as in the previous section. 

The Dispersion-Relation-Preserving (DRP) method developed by Tam and Webb 7 has 
been used to discretize these equations. The DRP scheme is an optimized fourth-order 
central finite difference scheme with a seven-point stencil in space and an optimized second- 
order multistep scheme in time. Non-reflecting boundary conditions are needed at the 
outer boundaries of the computational domain. The asymptotic non-reflecting boundary 
conditions of Tam and Webb 7 have been used for the two-dimensional computations in 
this paper and they have also been extended for the three-dimensional cases. The complete 
description and derivation of the schemes and the non-reflecting boundary conditions are 
given by Chung 1 . 


4. Acoustic Reflection and Scattering by Two-Dimensional Bodies 

In this section a number of numerical examples are given of acoustic reflection and 
scattering by 2-d bodies. Solutions are obtained using the IMM and are compared with 
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either analytical solutions or numerical solutions obtained using the solid wall boundary 
conditions proposed by Tam and Dong 6 . 

4.1. Acoustic Reflection of A Single Initial Acoustic Pulse By An Infinite Flat Solid 
Wall 



Fig. 4. Sketch of Computational Domain for Infinite Wall Reflection Problem. 

First, the acoustic reflection of a two-dimensional acoustic pulse by a plane wall is 
considered. This is shown schematically in figure 4. The wall is located at y = 0. This is 
an initial value problem. An acoustic pulse is generated by an initial pressure disturbance 
with a Gaussian spatial distribution. The initial conditions are: 

p = p = exp{-^[x 2 + (y- 25) 2 ]}, u = v = 0 (4.11) 

25 

The source is placed at (0,25). The fluid domain is 201 by 201. There is uniform mean flow 
with Mach number 0.5 parallel to the wall. The time step, At is 0.05. 

In order to simulate the infinite wall using the IMM, an extra wall region is needed as 
shown in figure 4. The thickness of this wall region is chosen to be 40. This thickness could 
be smaller, but in that case the source would be too close to the non-reflecting boundary 
of the total domain, and some wave reflections would occur at the bottom boundary. An 
extra wall region is needed only such wall situations at the boundary of the domain, not 
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for scattering by finite objects. The computations are then carried out directly. No stencil 
change is needed, no solid wall boundary conditions are implemented. The presence of the 
wall does not affect the speed of computations. Even though the extra wall region increases 
the amount of computations, the overall computing time is decreased compared to that 
using traditional solid wall boundary conditions. 
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Fig. 5. Pressure Contours at t — 100 for Reflection of an Acoustic Pulse by an Infinite Wall With Mean 
Flow Mach Number M = 0.5. 



Fig. 6. Pressure distribution on y axis and x = 50 at time t = 100 for Reflection of an Acoustic Pulse by an 
Infinite Wall With Mean Flow Mach Number M = 0.5. 

Figure 5 shows a calculated pressure contour associated with the acoustic pulse at t — 
100. At this time, the pulse has reached the wall and has been reflected off the wall creating 
a double pulse pattern. The entire pulse has been translated downstream by the mean flow. 
Figure 6 shows the corresponding computed pressure waveform along the line x = 50, which 
passes the center of the pulse. The numerical result is compared with the analytical solution 
and the agreement is good. But some small errors can still be seen in the reflected waveform, 
both in its amplitude and phase. The amplitude error is mainly due to the finite choice of 
mean density ratio. This was discussed in section 2. The phase error is caused by the fact 
that the location of the wall can not be defined exactly. This is discussed further in the 
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following example. 

4 . 2 . Acoustic Reflection of Time Periodic Acoustic Waves By An Infinite Flat Solid 
Wall 

The reflection of a periodic acoustic wave train by a solid wall is considered in the absence 
of a mean flow. The computational domain and the locations of the wall and the source 
are the same as in the previous example. The acoustic wave train is generated by a time 
periodic source in the energy equation. The source takes the form 

H = 0.01exp{^|p-[a: 2 + {y - 25) 2 ]} cos(wi) (4.12) 

u) is the angular frequency. 10 points per wavelength are used, so that u) = 7r/5. 


p 

0 001 

000077777 # 
0 0006U6M 
OOOOOMU* 
0000111111 
-0.000111111 


-0 000777776 
•0.001 



Fig. 7. Pressure Contours at t = 180 for Reflection of Acoustic Wave Train by Infinite Wall. 

The simulation is carried out with zero initial conditions. After the transient solution 
has propagated out of the computational domain, the pressure fluctuation is time periodic 
with an angular frequency u. Figure 7 shows the computed pressure contours adjacent to 
the solid wall at time t = 180 in the left half of the domain. The interference pattern is 
due to cancellation between the incident and the reflected waves. Figure 8(a) gives the 
corresponding pressure waveform along the y-axis. The analytical solution is also plotted. 
Noticeable errors can be seen in both the amplitude and the phase. In this case the wall is 
located at y = 0 in the analytical solution. Figure 8(b) also shows the pressure waveforms 
at the same time, but in this case the position of the wall is at y = 0.5 in the analytical 
solution. It can be seen that the agreement between the computed and analytical solutions is 
much better. This simulation demonstrates that there is a one grid spacing (Ay) inaccuracy 
in the definition of the wall position in the IMM. The mean density ratio in the IMM is 
specified in the following way: po — 1/30 when y < 0; po = 1 when y > 1, so the wall could 
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(a) (b) 

Fig. 8. Pressure Distribution on y axis at x = 0 and t = 180 for Reflection of Acoustic Wave Train by 
Infinite Wall. For Analytical Solutions, Wall is at (a) y = 0.0, (b) y = 0.5. 


be anywhere between y — 0 and y — 1 and the numerical solution could be unchanged. 
This is disadvantage of the IMM. If the source is not too close to the wall, or enough grid 
points are used between the source and the wall, then this error is in the acceptable range. 
But this disadvantage can have a positive effect when curved solid boundary problems are 
solved using the IMM. This is shown in the following example. 

4.3. Acoustic Scattering of Time Periodic Acoustic Waves By A Finite Flat Plate 
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Fig. 9. Sketch of Computational Domain for Finite Plate Scattering Problem 

The scattering of a periodic acoustic wave train by a thin flat plate of finite length is 
considered as shown in figure 9. The length of the plate, L — 25. The domain is 8 L by 8 L. 
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Fig. 10. Pressure Contours at t = 194 for Scattering of Periodic Source by a Finite Plate 



Fig. 11. Pressure Distribution Along Upper Boundary of Computational Domain for Scattering of Periodic 
Source by a Finite Plate 


The plate is located at the center of the computational domain. The acoustic waves axe 
generated by a simple harmonic source located at a distance of L from the middle point of 
the plate. The source term is incorporated into the energy equation and has the form 

H = 0.01 ex P{ — 777^2 fc 2 + (y ~ 25 ) 2 ]} cos(u>f) (4.13) 

( 1 -/ 0 ) 

In this numerical simulation, L = 4A, so u> = 2n/\ = 0.327T. This is only approximately 6 
points per wavelength. 

The computations are conducted using two methods: the IMM and the solid wall bound- 
ary conditions developed by Tam and Dong 6 . The thickness of the plate is Ay in the IMM. 
and zero in the solid wall boundary condition method. The pressure on the two sides of the 
plate is different, so in the IMM at least one Ay of thickness is needed. That is, in two rows 
of length L, the mean density equals 1/30. When an acoustic wave train impinges on the 
plate, the wave is scattered. In the shadow region behind the plate, acoustic waves radiated 
directly from the source are blocked so the sound pressure consists only of the contributions 
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from waves diffracted by the two sharp ends of the plate. Figure 10 shows the pressure 
contours computed using the IMM at time t = 194. The diffraction pattern behind the 
plate and the scattering pattern in front of the plate can be seen clearly. Figure 11 shows 
the corresponding pressure distributions along the upper boundary of the computational 
domain obtained from both the IMM and Tam and Dong’s method 6 . The agreement be- 
tween the two solutions is very good. The computing time for the IMM is two-thirds that 
for Tam and Dong’s method. Also, the coding in the IMM is extremely simple. In the 
IMM, the amount of coding work and the computing time do not change at all when there 
is an solid object present. This is one of the major advantages of the IMM. 

4.4. Acoustic Scattering of Time Periodic Acoustic Waves By A Circular Cylinder 



Fig. 12. Sketch of Computational Domain for Circular Cylinder Scattering Problem 

The scattering of a periodic acoustic wave train by a circular cylinder is considered next. 
This is sketched in figure 12. The computational domain is 201 by 201. The cylinder is 
placed at the center of the domain and has a diameter D = 25. The acoustic wave train is 
still generated by a time periodic source in the energy equation. The source and its location 
are the same as those described in section 4.2. 

The important difficulty in this example is how to deal with the curved solid bound- 
aries. Uniform Cartesian grids and high-order finite difference DRP schemes have been 
used throughout this paper. The use of a uniform grid has advantages in the maintenance 
of good dispersion and dissipation properties, as does the high-order finite difference DRP 
scheme. However, their use also presents a problem. This is the non-conformity of the grids 
with the boundaries of curved bodies. Kurbatskii and Tam 4 have recently introduced a 
technique for Cartesian treatment of curved walls for high-order finite-difference schemes. 
In the present calculations we used the nearest grid points to the boundary to define the 
body shape. So the curved boundary is approximated by a staircase boundary. This results 
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Fig. 13. Pressure Contours at t — 180 for Scattering of Periodic Source by a Circular Cylinder. 

D = 25. 



X 

Fig. 14. Pressure Distribution Along Upper Boundary of Computational Domain for Scattering of Periodic 
Source by a Circular Cylinder 


in some errors in the computation. But, as demonstrated in the following example, this 
error may be small if the IMM is used. For the circular cylinder, the wall region is defined 
as ( x 2 4- y 2 ) 1//2 < D/ 2. Inside this region, the mean density is set equal to 1/30. This is all 
that is needed in order to define the presence of the cylinder. After this the computations 
are carried out directly. It is obvious that the coding work and computing time are the 
same with or without the object. Also, there is no difference at all if the boundary is curved 
or flat. This is the most important advantage of the IMM. 

Figure 13 shows the pressure contours computed at time t = 180. The scattering pattern 
behind and in front of the cylinder is seen clearly. Figure 14 shows the pressure distribution 
along the upper boundary of the computational domain. The analytical solution is also 
shown. The agreement between the computed and analytical results is good. It can be 
concluded that the use of a stair-stepped surface to approximate the curved boundary is a 
reasonable approximation in the IMM. Since in the IMM the position of the object surface 
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could be anywhere in a range of one grid spacing, this has a smoothing effect on the staircase 
surface and makes the surface approximation more accurate. It can be seen that the IMM 
is a very efficient and convenient method to treat curved boundaries. 

5. Acoustic Reflection and Scattering From 3-D Bodies 

In this section, two problems of acoustic scattering by three-dimensional bodies are solved 
using the IMM. The first problem is the acoustic scattering of a time periodic acoustic wave 
by a sphere. This is an axisymmetric 3-d problem. The analytic solution is given by 
Morris 5 and is compared with the numerical solutions. The second problem is acoustic 
radiation and scattering from a cylindrical shell with a acoustic source placed inside the 
shell. No analytical solutions are available for this problem, so that only numerical results 
are presented. 

5.1. Acoustic Scattering of Time Periodic Acoustic Waves By A Sphere 



Fig. 15. Sketch of Computational Domain for Sphere Scattering Problem 

The scattering of a periodic acoustic wave train by a sphere is considered as shown in 
figure 15. The domain is 51x51x51. The sphere is placed at the center of the computational 
domain. This is also the origin of the coordinates. The sphere diameter is D = 12.5. The 
acoustic wave train is generated by a time periodic source in the energy equation. The 
source term has the form 

H = 0.01 ex P{ ( 2 ^ [* 2 + (y + 16 ) 2 + z 2 ]} 

cos (u;t) 


(5.14) 
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u) is the angular frequency. 10 points per wavelength are used, so that u) = 7r/5. The center 
of the source is at (0,-16,0). 



Fig. 16. Pressure Contours at t = 45 for Scattering of Periodic Source by a Sphere. 

The sphere is the wall region, which is defined by \Jx 2 + y 2 + z 2 < D/2. Inside this 
region, p 0 — 1 /30. It is clear that the boundary of the sphere is approximated by a staircase 
boundary. Figure 16 shows the computed pressure contours at the z = 0 section at time 
t = 45. This section passes through the center of the sphere. The scattering pattern behind 
and in front of the sphere is clear. Figure 17 is the pressure waveform along the x-axis at 
this section and at the same time. Figure 18 is the pressure waveform along the j/-axis. 
The numerical results may be compared with the analytical results and good agreements 
are achieved. Most of the disagreement occurs because of the relatively coarse grid used in 
this simulation. 

This numerical example demonstrates that the IMM is applicable to 3-d acoustic scatter- 
ing problems. The coding work would be tremendous if the traditional solid wall boundary 
conditions were used for 3-d scattering cases. For the IMM, there is no extra coding work 
and computation time at all when the 3-d scatterers exist. 

5 . 2 . Acoustic Radiation and Scattering of Time Periodic Acoustic Waves Inside A 
Cylindrical Shell 

As a final example, the radiation and scattering of a periodic acoustic wave train inside 
a cylindrical shell is considered. The acoustic wave train is generated by a time periodic 
source in the energy equation, the same as that in section 5.1, located at (0,0,0). The 
domain is 51x51x51. The shell is placed at the center of the domain. The inner radius of 
the shell is r\ = 6.25, the outer radius is r 2 = 9.25. The thickness of the shell is 3. The 
length of the shell is h = 26 and the axis of the shell aligns with z-axis. The computational 
domain is sketched in figure 19. 

The wall of the shell is defined as rq < \Jx 2 + y 2 < r 2 , and \z\ < 13. Inside this region 
p 0 = 1 /30. So the inner and outer boundaries of the shell are all approximated by staircase 
boundaries. As an example of the calculated pressure contours figure 21 shows the pressure 
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Fig. 17. Pressure Distribution Along x-axis at time t = 45 for Scattering of a Periodic Source by a 

Sphere. 



Fig. 18. Pressure Distribution Along y-axis at time t = 45 for Scattering of a Periodic Source by a 

Sphere. 


contours at section y = 0, which cuts through the center of the shell. It can be seen that 
when an acoustic source is placed inside the shell, the acoustic waves can radiate from the 
two open ends. These two open ends diffract waves and act like sources for the external 
acoustic field. Outside the shell the pressure field is stronger closer to the ends, and weaker 
away the ends. Inside the shell a standing wave pattern can be seen 

This numerical simulation demonstrated that the IMM can deal with various kinds of 
3-d objects very easily. The complex geometries, so long as they reasonably smooth, do not 
represent any difficulty in the IMM. 

6. Conclusions 

In this paper we have introduced an efficient method to simulate solid wall boundaries, 
the Impedance Mismatch Method. This method was applied to a number of 2-d and 3-d 
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Fig. 19. Sketch of Computational Domain for Shell Scattering Problem 


scattering problems with both flat and curved boundaries. Many advantages of the IMM 
have been found and demonstrated. No special solid wall boundary conditions need to 
be implemented. No stencil changes are involved because of the presence of solid objects 
and the coding is very easy. The computations are much faster than when the traditional 
solid wall boundary treatments are used. There is no difficulty for any reasonable smooth 
geometry. No matter whether the solid boundary is fiat or curved, the amount of coding 
work and computing time are the same. Some disadvantages of the IMM have also been 
revealed. The accuracy of the computations depends on the value of the mean density ratio 
in the wall region. Also there is an one-grid-size error in the wall position. An extra wall 
region is needed for the infinite wall case. Finally, a staircase approximation is used to 
approximate curved boundaries. This paper has demonstrated that the IMM is a promising 
method for the simulation of acoustic reflection, scattering, and diffraction problems. 
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Abstract 

Supersonic rectangular jet flow and far field noise pre- 
dictions are made by solving the governing equations 
using advanced numerical techniques on parallel pro- 
cessors. The computational domain begins at the jet 
nozzle exit and contains the jet plume and a small re- 
gion of the acoustic near field. The equations solved 
for the interior grid points are the full 3D Navier Stokes 
equations. The far field boundary points are determined 
by unsteady, nonlinear characteristic based nonreflect- 
ing conditions. To model the jet nozzle exit flow, a set 
of equations axe developed to simulate many features 
of this flow that have been experimentally observed to 
influence the jet and its radiated noise. A Kirchhoff 
method is used to determine the far field noise from in- 
formation extracted from the finite computational do- 
main. Each set of governing equations is spatially dis- 
cretized by a sixth order central difference scheme and 
advanced in time using fourth order Runge-Kutta in- 
tegration. Spurious high wave number fluctuations are 
damped by a nonlinear dissipation algorithm that has 
a minimal effect on the acoustic solution. The code has 
been efficiently implemented on the CMS using CMFor- 
tran (essentially HPF) and should be easily ported to 
platforms running HPF (such as the SP2). Numerical 
results indicate that the algorithm, which contains no 
model constants (aside from the nozzle exit conditions), 
is capable of reproducing many experimentally observed 
rectangular jet flow and noise features. 

1 Introduction 

Jet noise analysis and reduction have been topics of 
research since the introduction of the jet engine as 
a propulsion device for aircraft during World War II. 
Many tools have been developed and incorporated into 
the design process to reduce the annoyance of jet noise. 
The Federal Aviation Administration began placing 
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strict regulations on the noise produced by aircraft 1 and 
this has caused a renewed interest in noise prediction in 
the scientific community. This is particularly true in 
light of the national interest to develop the High Speed 
Civil Transport (HSCT). Noise reduction is considered 
to be a crucial technology required for a viable design. 

The prediction strategy that is currently receiving the 
most attention is direct simulation. Jet noise research 
has included many different implementations of this ap- 
proach. 2-5 The implementation under investigation in 
this paper solves the 3D full Navier Stokes equations in 
a domain that includes the noise source region and a 
small portion of the acoustic field. The acoustic near 
field solution is then used as an input to a Kirchhoff 
method to determine the fax field noise. 

Consistent with experimental evidence, 6,7 the noise 
source is assumed to be dominated by the evolution of 
large scale coherent structures in the jet shear layer and 
thus only these scales are resolved. The influence of the 
small scales is assumed to be represented by numerical 
dissipation. Thus, no turbulence model is employed and 
consequently, there are no adjustable model constants 
used in the interior domain. 

The computational domain in which the Navier 
Stokes equations axe solved begins at the jet nozzle 
exit plane. A model is therefore required for the nozzle 
exit flow and can have an appreciable effect on the nu- 
merical solution. The effects of nozzle exit conditions 
on experimental and numerically simulated jets have 
been studied by many investigators. 8-15 Hussain and 
Husain 8 found experimentally that the development of 
the jet depends on the nozzle boundary layer momen- 
tum thickness distribution. The azimuthal variation has 
been shown by them to produce noticeable effects on the 
spreading rate of elliptic jets. These effects are the re- 
sult of the influence of the momentum thickness on the 
generation of coherent structures. King et al. 12 found 
that nozzle imperfections as small as 0.2% of the nozzle 
exit diameter may have a significant effect on the de- 
velopment of supersonic axisymmetric jets. They used 
this information to develop methods of enhancing jet 
mixing. The initial turbulence intensity was found have 
a significant effect on the turbulence amplification rate 
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in the near-field region of the jet by Grinstein et al. 11 
Quinn 14 found differences in the spreading rate of two 
jets operating under essentially the same flow conditions 
and geometry. These discrepancies were attributed to 
facility differences. 

This experimental evidence suggests that the salient 
features of a laboratory facility will have an apprecia- 
ble effect on the development of a jet. These effects 
can be observed in the form of varying potential core 
lengths, turbulence levels and jet spreading rates. To 
compound this potential problem, the nozzle exit con- 
ditions of the rectangular jet simulated in this paper 
are not defined precisely. The measurement techniques 
available to Kinzie 16 did not permit a comprehensive 
study of the nozzle exit. Due to this uncertainty, a gen- 
eral model for the nozzle exit flow has been developed 
that can turn on or off some of the features that have 
been observed to significantly influence the jet develop- 
ment. The results presented in this paper are confined 
to studying the effect of modal excitation. The effects 
of nozzle exit turbulence levels and corner vortices are 
discussed in Chyczewski(1996). 17 

Since the computational domain is limited to just a 
small region of the acoustic near field, a method is re- 
quired to extrapolate the solution to the far field. In 
this work the Kirchhoff method is employed. 18 It con- 
sists of constructing a surface S on which the acoustic 
solution can be reliably calculated. The acoustic solu- 
tion at any location outside of this surface can then be 
determined by the Kirchhoff formula. The solution is 
exact for sound radiation outside of a surface S if that 
radiation is governed by the convective wave equation. 
However, in the application of this method to the jet 
noise problem, finding such a surface is difficult. This 
issue has been addressed previously by some investiga- 
tors who have found that the far field solution is not 
very sensitive to the location of the surface if some pre- 
cautions are taken. Lyrintzis and Mankbadi 19 found 
that placing the surface at least one diameter away from 
the jet centerline is sufficient to obtain accurate solu- 
tions. Freund et al. 20 performed a study analyzing the 
effects of using open Kirchhoff surfaces and found that 
it does not introduce significant errors if the surface 
passes through the region between the noise source and 
the observer. 

In the next section, the governing equations axe de- 
scribed. This consists of discussing the specific form of 
the Navier Stokes equations, presenting the boundary 
equations, which includes the model nozzle exit condi- 
tions, and finally presenting the Kirchhoff formulation 
used here. In section 3 the numerical approach is out- 
lined. Special attention is given to the artificial dissi- 
pation model. Next, in section 4, rectangular jet noise 
prediction results are presented and compared to exper- 
imental data. Finally, in section 5, some conclusions are 
drawn. 


2 Governing Equations 

A supersonic rectangular jet flow is a nonlinear, vis- 
cous, unsteady, 3D problem. As such, it is governed by 
the full, compressible, 3D Navier Stokes equations. A 
nondimensional conservative form of these equations is 
used in this work (see Hoffmann(1989) 21 ). By them- 
selves, these equations are not sufficient to model the 
jet problem. Boundary conditions are required to al- 
low flow and acoustic waves to pass through the far 
field boundaries of the computational domain as well as 
to model the flow entering the domain from the nozzle 
exit. This section presents these equations as well as the 
Kirchhoff formulation used to extrapolate the acoustic 
solution to the far field. 

2.1 Nonreflecting Boundary Conditions 

Several approaches to the specification of nonreflect- 
ing conditions at far field boundaries have been devel- 
oped. These approaches can be classified into three cat- 
egories: asymptotic solutions, 22,23 Fourier decomposi- 
tion 24,25 and quasi one-dimensional analysis. 26,27 The 
most recent set of conditions based on the asymptotic 
solution of the linearized Euler equations are due to 
Tam and Webb. 23 These asymptotic conditions have 
been quite successful at reducing boundary reflections 
for many model problems. 

Giles 24 derived approximate unsteady boundary con- 
ditions for two-dimensional problems by performing a 
Fourier decomposition of the linearized Euler equations. 
They have been applied to turbomachinery problems by 
Giles and have been found to be effective. When imple- 
mented with a buffer zone, these conditions have been 
able to permit nonlinear vortical structures to leave a 
computational domain with little reflection. 25 

A drawback of both the asymptotic and Fourier meth- 
ods is that their derivation employs a set of equations 
that have been linearized with respect to a reference so- 
lution. In many cases, such as the rectangular jet prob- 
lem under consideration here, the reference solution is 
not known a prion and must be developed as the equa- 
tions are integrated. Experimentation with the rectan- 
gular jet problem suggests that asymptotic and Fourier 
methods are not capable of establishing a reasonable 
reference, or time averaged, solution when the initial 
condition is a quiescent fluid. 

Given this difficulty, the quasi one-dimensional 
boundary procedure developed by Thompson 26,27 is 
employed. The approach consists of decomposing the 
full nonlinear Euler equations into modes of definite ve- 
locity and specifying nonreflecting conditions for those 
modes that have a velocity directed into the computa- 
tional domain. These conditions have been shown to be 
able to allow large amplitude disturbances to leave the 
domain with little reflection. 26 
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2.2 Nozzle Exit Conditions 

A complete prescription of the nozzle exit conditions 
requires the specification of both steady and unsteady 
characteristics. They are described in the following two 
subsections. 

2.2.1 Steady Nozzle Conditions 

A nearly uniform velocity profile was found at the noz- 
zle exit by Kinzie. 16 This indicates that viscous effects 
are confined to locations very close to the nozzle wall. 
Thus the simulated jets use uniform profiles for density, 
axial velocity and pressure. Accurate measurements of 
momentum thickness variations that may exist around 
the nozzle lip were not performed and are thus not ac- 
counted for in the nozzle model. 

The values of the exit variables are found from the 
experiment. The exit Mach number, Afj, is 1.54, the 
acoustic speed of the jet is Cj = 0.82cqo and since the 
jet is ideally expanded, the jet exit pressure is the same 
as the ambient pressure {pj — Poo)- From this informa- 
tion, the steady exit density and velocity can be found. 
In this paper, the nozzle exit flow is assumed to be 
purely axial. The effects of lateral exit flow compo- 
nents induced by nozzle exit corner are considered in 
Chyczewski(1996). 17 

2.2.2 Unsteady Nozzle Conditions 

There is a very limited amount of information available 
in the literature that discusses the unsteady features of a 
supersonic nozzle exit flow. In fact, the authors have not 
seen any published data that characterizes the unsteady 
features to the extent that is required to reproduce the 
nozzle exit conditions completely. This is most likely 
due to the extreme difficulty of collecting such data. 

To compensate for this lack of information, a general 
set of unsteady nozzle exit conditions have been devel- 
oped that can specify the disturbance spatial distribu- 
tion, amplitude, temporal behavior and phase relation 
around the nozzle lip. By controlling the phase relation, 
different modes (flapping or varicose) can be excited at 
the nozzle exit. This is similar to the artificial excita- 
tion used by many experimentalists. 7,16,28-30 With this 
model, many different features can be investigated. In 
this paper, the investigation is confined to studying the 
modal excitation. 

The velocity perturbations are calculated from the 
following relation : 

4 2 

u\ v * or ^ = + ^ + ft) 

t=i i=i 

The contributions of each of these terms is given in the 
following sections. 


Temporal Behavior. 

The inner summation in equation 1 is over contri- 
butions from two characteristic frequencies. These two 
frequencies are the screech tone frequencies found in the 
minor axis plane of the experimental jet. 16 These fre- 
quencies can also be determined using the linear shock 
ceil model and weakest link theory developed by Tam. 31 
For our problem these frequencies are /i = 9606 Hz 
(Strouhal number = 0.31) and p 2 = 26367 Hz (Strouhal 
number = 0.84). 

A random component to the excitation is supplied by 
0” in equation 1. It is initialized to zero at the begin- 
ning of each run. It is then updated at each timestep, 
n, by the following : 

#* = #*■‘±0 ( 2 ) 

The amplitude of the phase shift between time steps 6 
is 5.4 degrees. This value was found to give broad fre- 
quency spectra with an upper band limit that is near 
the highest frequency that the grid and scheme can re- 
solve. 

Spatial Distribution. 

The perturbation on the entire nozzle lip is deter- 
mined by building it up from the contributions of the 
four walls. This is done in equation 1 by the outer sum- 
mation. is the spatial amplitude function for each of 
the walls. It is a Gaussian function centered on the lip 
line of each wall. The half width of the Gaussian is one 
tenth of the short dimension of the nozzle. The function 
is tapered to zero amplitude near the corners of the noz- 
zle. This spatial function is selected since a Gaussian 
is a representative distribution for wall bounded shear 
layer perturbations (see Kinzie(1995) 16 for example). 

Mode Excitation. 

Two different modal excitations are considered in this 
paper. These axe the varicose and flapping modes and 
have been found in the experimental jet by Kinzie. 16 
The varicose mode is characterized by symmetric shed- 
ding of coherent structures from the nozzle lip. In con- 
trast, the vortex shedding is asymmetric for the flapping 
jet case. This is illustrated in figure 1. These different 
modes are excited in the jet by controlling phase dif- 
ferences between the walls of the nozzle. The phase 
difference is controlled by the angle in equation 1. 
Whether or not there is a velocity component contribu- 
tion to the perturbation from a wall is determined by 
the parameter Cj. The values of these two parameters 
^ for each of the velocity components is given in tables 1 
y and 2 for the varicose and flapping modes. 

The final parameter in equation 1 is a. It speci- 
fies the peak amplitude of the velocity perturbation. A 
value of 0.02 is used and corresponds to an RMS fluc- 
tuation level near 1.36 percent of the exit velocity. The 
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Figure 1: Vortex shedding in varicose (top) and flapping 
(bottom) excited jets. 
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Table 1: Values of c* and 0 t for varicose mode excita- 
tion. 


RMS levels vary slightly from run to run due to the 
random nature of the excitation. 

The instantaneous velocities are obtained by adding 
the steady contribution to the perturbation. Given 
these velocities, the pressure and density are found from 
conditions of constant total temperature and entropy. 
The first condition was found to be reasonable for a jet 
flow under these conditions by Troutt and McLaugh- 
lin. 32 The assumption of isentropic excitation is jus- 
tified since these perturbations most likely originated 
from acoustic disturbances upstream of the nozzle exit. 


2.3 Kirchhoff Formulation 


The moving surface formulation given by Farassat and 
Myers 18 for a rigid surface in rectilinear motion is em- 
ployed. It gives the acoustic pressure p l at location x 
and time t as a function of the pressure on a suitably 
defined surface 5 (where x is outside of S ): 


4.p'(f,() = // s [ RT |! S3 

+ [[\ ill 

JJs [r 2 (l — M r 


dS 

dS 


where 


r = \f\, f=x — y(r), M r = M • r/r, 


(3) 

(4) 

(5) 


Er = -n ■ Vp' + (M ■ n)(M ■ Vp') 
cos 8 - M ■ n M • ft dp f 
+ Cecil -M r ) Coo dr ’ 
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Table 2: Values of C{ and 0 l for flapping mode excita- 
tion. 


52 = (l-M r ) 2 ^ COS0 -M-n) (8) 

M is the Mach number of the moving surface which 
for the static jet problem is zero. The vector r is the 
vector difference between the observer location and the 
location of the Kirchhoff surface element (it varies with 
each location on the surface), n is the normal vector 
pointing out of the Kirchhoff surface, the angle 8 is 
measured between the vectors r and n, and c Q 0 is the 
freestream sound speed. The integrands are evaluated 
at the Kirchhoff surface emission time r * which, for a 
stationary surface, is given by 

r* = t - r/coo (9) 

A complete description of the coupling of the Kirchhoff 
method into the Navier Stokes code is given in Ozyoruk 
and Long 33 and is therefore not described here. 

3 Numerical Algorithm 

The governing equations are discretized in a finite differ- 
ence context using fourth order accurate Runge-Kutta 
time integration and sixth order accurate spatial dis- 
cretization. The computational domain is illustrated in 
figure 2 which also shows the coordinate system. The 
center of the nozzle exit is located at (x, y, z) = (0, 0, 0). 
Details on the grid generation strategy can be found in 
Chyczewski and Long(1995) 34 and Chyczewski(1996). 17 


3,1 Artificial Dissipation 

A desirable feature of this numerical algorithm is the 
explicit control of the amount of dissipation applied to 
the scheme. Unlike upwind methods, Runge Kutta - 
central difference techniques contain very little implicit 
dissipation. Instead, explicit filters are used for stabil- 
ity and to prevent odd-even decoupling errors. Selec- 
tion of an appropriate dissipation scheme is paramount 
in calculations where one wishes to extract the low am- 
plitudes and high frequencies associated with acoustics. 
Jameson et al. 35 proposed a blend of split second and 
fourth order dissipation. While this dissipation is ro- 
bust near discontinuities, it may significantly contami- 
nate the acoustic solution. 
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Figure 2: Three dimensional view of the grid. 

A combination of second and sixth order dissipation is 
used here. This dissipation has been applied to the non- 
linear acoustic benchmark problems 34 and has proven 
capable of propagating acoustic waves in the presence 
of strong discontinuities. It is applied to the scheme as 
a correction to the residual and takes the following form 

T>(Q ) itjtk = T>t(Q)ij y k + IVQkM + ®c(Q)u.* ( 10 ) 

Each term is determined as follows (using as an ex- 
ample) : 


jet exit. Since the static pressure is essentially uniform 
there (ideally expanded jet), using the static pressure is 
inappropriate. 

The second order dissipation coefficient, is ex- 
plicitly set to zero in regions where the flow is suffi- 
ciently smooth, i.e. v is below a specified value. To 
minimize the size of the region with second order dissi- 
pation, this value should be raised to a maximum that 
results in a stable scheme. An appropriate value has 
been determined to be 0.0050 by numerical experimen- 
tation. The locations where the second order dissipation 
is necessary is found to be confined to certain locations 
in the jet core and should have no effect on the radiated 
acoustic solution. This is illustrated in figure 3 where 
a snapshot of the locations in the minor plane where 
second order dissipation is applied at a typical instant. 



Z><(Q )>,M = e (2) 2?|(Q )ij,k + e (6) 2>f(QkM (11) 


where 


Figure 3: Snapshot of the locations in the minor axis 
plane where second order dissipation is applied. 


= (Q<+1 Ji* ~ Qi-l,j,k) (12) 

and 

(Q )i,j,k = ( — 20Qi j,* + 15(Qt+i,j,* + Q 

+ Q*-2j,Jfc) + (Qi+3,j,fc + Qi— 3,j,fc)) (13) 

The coefficients are determined in a manner very similar 
to that used by Jameson: 


t < 2) = fc (2) max(i/ j _i,i/ i ,i/ i+ i) 


(14) 


and 


e( 6) = max(0,fc <6) -e| 2) ) (15) 

where the values of and k < 4) used here are 1.5/4 
and 1.5/256, respectively. The flow gradient sensor v is 
given by 

{ 0i — 1 , j , k — 2fft ,j,k ~f~ 0i+l j,fc| 


Vi = 


Pi-1, j,k + 2 0i,j,k + Pi+l,j,k 


(16) 


0 is set equal to the total pressure. The Jameson 
scheme 35 uses the static pressure. For our jet calcu- 
lations, the highest flow gradients are found near the 


4 Results 

In this section, rectangular jet simulation results are 
presented and compared to experimental data. The re- 
sults are obtained by executing a run that consists of 
three phases. Since the initial condition of the sim- 
ulation is a quiescent fluid, one phase is required to 
allow transients to leave the domain and establish the 
jet. When the domain is free of transients, the Kirch- 
hoff integration is started. There is a transient conver- 
gence period required by Kirchhoff methods which is 
based on the furthest distance between the Kirchhoff 
surface and an observer location. 33 The time required 
to converge the Kirchhoff solution constitutes the sec- 
ond phase. The final phase is used to sample variables 
in the jet and in the near and far acoustic fields. Only 
the data collected in the last phase has been used to 
generate the results presented here. 

Figure 4 shows the centerline velocity distributions 
for the varicose and flapping excited jets and compares 
them to the experimental data. 16 The varicose jet sim- 
ulation results are shifted -2 D eq to match the potential 
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core length of the experiment. This common proce- 
dure is used so that the decay rates of the centerline 
velocity can be compared directly. This shift is applied 
to all comparisons of the varicose excitation simulation 
with experimental data presented throughout this pa- 
per. The flapping excited jet does not require a shift. 



Figure 4: Centerline velocity distributions. 


The comparison between the simulations and experi- 
ment after the end of the potential core shows that the 
simulation overpredicts the turbulent mixing slightly for 
both of the mode excitations since their centerline ve- 
locities decay at a faster rate. A distinguishing feature 
between the two simulation profiles is that the flapping 
excitation case has noticeable oscillations in the poten- 
tial core region. These oscillations, also present in the 
experimental data, are due to a shock cell structure. 
They are also present to a lesser degree in the vari- 
cose excitation case profile. In the experimental jet, 
this structure most likely originates in the throat of the 
nozzle as a result of an imperfect nozzle design. 

Recall that the static pressure prescribed by the 
steady nozzle exit conditions is set to the ambient pres- 
sure; however, superimposed on this steady condition 
are perturbations. These perturbations are likely re- 
sponsible for the weak shock cell structure found in the 
simulated jet. Why the flapping mode excitation pro- 
duces a stronger shock cell structure is not understood. 
It is interesting to note, however, that the simulation 
reproduces fairly accurately the amplitude, wavelength 
and phase of the shock cell structure (given the limited 
resolution of the experimental data). This may suggest 
that the mechanisms producing the shock cell struc- 
tures in both jets are similar and that the experimental 
shock cell structure is not solely due to an imperfect 
nozzle design. 

Since the evolution of large scale turbulent structures 
plays such an important role in supersonic jet noise gen- 
eration, some of the simulated properties of these struc- 
tures have been determined and compared to experi- 
mental data. In his experiment, Kinzie 16 used hot wire 
anemometry to measure the fluctuations in the shear 
layer. In a compressible flow, these wires are sensitive 


to the mass flux that is normal to the wire. Measure- 
ments were made in both the major and minor axis 
planes of the jet. When performing major axis plane 
measurements, the hot wire was oriented parallel to the 
minor axis plane, i.e., it was parallel to the wail from 
which the shear layer was emanating. This was done to 
improve the resolution of the shear layer measurements. 
When arranged in this manner, the normal component 
of the mass flux consists of two velocity components, 
i.e., the hot wire is sensitive to : 

m = y/(pu) 2 + ( pw ) 2 (17) 

Figure 5 shows the axial development of the RMS of 
the variable m defined in equation 17 (normalized by 
the jet exit mass flux) in the major axis plane. The val- 
ues plotted in the figure are the maximum RMS values 
through the shear layer for a given axial location. 



Figure 5: Maximum RMS levels of the mass flux normal 
to the hot wire in the major axis plane. 

The most apparent observation made from this fig- 
ure is that the simulation over predicts the peak am- 
plitude of the perturbations by approximately a factor 
of 2. Thus, there is significantly more turbulent mix- 
ing in the simulated jet compared to the experimental 
one. This is consistent with the centerline velocity decay 
profile presented earlier. A possible explanation for this 
discrepancy is the absence of a sub-grid-scale (SGS) tur- 
bulence model. In turbulent flow, the large scales of the 
turbulence are continually acted on by the finer scales. 
These fine scales behave as a dissipation mechanism for 
the larger scales. The algorithm used in this research 
does not address this issue explicitly. The algorithm 
applied here relies on the dissipation supplied by the 
numerical scheme to behave like the sub-grid scales of 
the true flow. The difference between the artificial and 
SGS dissipation is quantified in Chyczewski(1996) 17 and 
will not be discussed in detail here. It will just be men- 
tioned that a comparison of these two dissipation terms 
reveals that the SGS dissipation is usually larger than 
the artificial dissipation but the difference is not consid- 
ered significant enough to account for the discrepancies 
found here. The high amplitude perturbation found in 
the simulation may also be explained by the nozzle exit 
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conditions. Berman et al. 36 found that changing the 
exit conditions for their subsonic calculations can re- 
duce the peak amplitude of the perturbations. 

Aside from the peak amplitude discrepancy, a com- 
parison of the trends in figure 5 between the simulation 
and the experiment is encouraging. There is a high 
amplitude growth rate prior to the end of the poten- 
tial core, the value of which compares well between ex- 
periment and simulation. After this region, there is a 
saturation period and finally a gradual decay in the am- 
plitude. 

The power spectral density of the time series used to 
determine the RMS values discussed above for the vari- 
cose excitation in the major axis is presented in figure 
6. The sample used to determine these spectra spans 
the final phase of the run described at the beginning of 
this section. This phase consists of 32,768 timesteps. 
In order to make handling the data less cumbersome, 
a sample was taken once every 16 timesteps. This pro- 
cedure compromises no information since the sampling 
frequency is still much higher than the frequencies ex- 
pected to be produced by the jet. Thus, the length of 
the entire sample is 2048 steps. To reduce the errors as- 
sociated with using a finite (and relatively small) sample 
record, the 2048 sample is divided into 15 records that 
contain 256 steps. The intervals overlap one another by 
128 elements. For example, the first interval contains 
samples 1 through 256, the second contains 129 through 
384, and so on. The spectra presented in figure 6 are 
obtained by averaging the spectra obtained from each of 
the 15 intervals. A Hanning window is also used to re- 
duce the errors associated with the finite record length 
(see Bendat(1986) 37 for example). Since the area un- 
der the power spectral density is equal to the RMS of 
the fluctuation, which have been previously discussed, 
the experimental spectra have been scaled so that they 
have RMS values equal to the simulation. This allows a 
direct comparison of the spectral distribution of energy. 
The spectra are presented as a function of the Strouhal 
number, which is the ratio of the frequency to a char- 
acteristic frequency. The characteristic frequency used 
in these spectra is 31387 Hz, which is the experimental 
jet exit velocity divided by the equivalent diameter of 
the nozzle exit. 16 

The comparison between the experimental and simu- 
lation spectra is favorable. In both the experiment and 
the simulation the energy moves to lower frequencies at 
the larger axial locations. This also agrees with obser- 
vations made by Troutt. 29 The spikes in the major axis 
experimental spectra near a Strouhal number of 0.5 are 
due to the shock cell structure found in the jet. They 
are the screech tones due to the phase locking of the 
radiated noise from the interaction of the large scale 
structures with the shock cell structure and the excita- 
tion of instability waves at the nozzle lip. These tones 
are most likely absent from the simulated jet due to its 
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Figure 6: Power spectral density (PSD) of the mass flux 
in the at the maximum RMS location through the shear 
layer in the major axis plane for the varicose excitation. 
The numbers to the right of each spectra denotes the 
axial location. The heavy lines are the experimental 
results. 



Figure 7: Near field sound pressure level (SPL) contours 
in the minor axis plane of the simulated jet excited by 
a varicose mode. 

relatively weak shock cell structure. 

In figure 7 sound pressure level (SPL) contours are 
shown in the minor axis plane of the simulated jet ex- 
cited by a varicose mode. The definition of an SPL used 
here, consistent with Kinzie, 16 is : 

SPL = 201og 10 [§^] (18) 

*i ref 

where 

Pref = X 10 - 6 ) ^ ( 19 ) 

ifltm m 

The pressure in the anechoic chamber, P c h , during his 
experiments was 3080 N/m 2 . The atmospheric pres- 
sure, Patm, is 1.01325 x 10 5 N/m 2 . The root mean 
square of the pressure, P rma , is determined in the same 
manner as the data used for the hot wire comparisons. 
Consistent with the experimental data, 16 the noise ap- 
pears to be generated at the end of the potential core 
and is directed in the downstream direction. 

The remainder of this section will present and discuss 
far field noise predictions. For the far field noise predic- 
tions presented here the Kirchhoff surface is placed as 
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close to the jet as possible without intersecting regions 
of the jet plume where there are significant hydrody- 
namic fluctuations. Close proximity to the jet is desir- 
able since the grid resolution is the finest there. Also, 
the closer the surface is to the jet, the less the pres- 
sure waves have to travel before reaching the surface. 
Thus, the numerical damping and dispersion errors are 
minimized. It is important not to place the surface too 
close to the jet since the Kirchhoff formula may inter- 
pret some hydrodynamic fluctuations as radiating ones. 

The entire Kirchhoff surface is defined by four sur- 
face segments. Two surfaces are in the (i,j) plane and 
two are in the (t,fc) plane ( i is the axial index). The 
two end planes (perpendicular to the jet axis) are omit- 
ted. Since the noise propagation is predominantly in 
the downstream direction, the downstream end plane 
should not make any contributions to the noise levels 
at the observer locations used here (described below). 

The grid used in the present study (shown in figure 
2) is highly clustered near the nozzle lip region. This 
stretching varies in the axial direction so that a con- 
stant j or k grid line gradually moves away from the jet 
axis. This is very convenient in terms of deciding how 
to define the Kirchhoff surface. Each of the four surface 
segments can be placed very near the nozzle lip line at 
the nozzle exit plane. As the jet develops in the ax- 
ial direction, and the mixing increases, the surface will 
gradually move away from the jet axis so that it never 
does intersect high mixing regions. 

The Kirchhoff surface is illustrated in figure 8. The 
extent of the surface in the axial direction is from x — 
1 D eq to 20 D eq . The surface is not extended to the wall 
so that the possibility of any effects it may have on the 
acoustic solution is reduced. The surface terminates at 
20 diameters to reduce the effects of any reflections that 
may occur from the downstream boundary. The range 
of the domain covers the axial region where turbulent 
mixing noise is known to be generated, i.e., the region 
near the end of the potential core. 

The lateral locations of the surface were determined 
by considering the mass flux perturbation levels. Fig- 
ure 9 shows contours of the RMS of such perturbations 
in the major axis plane of the jet excited by a varicose 
mode. The grid (showing every tenth i grid line) is over- 
laid on the figures. The arrows (between x = llD eq and 
12 D eq ) point to the grid lines that seem to be optimal 
for the planes of the Kirchhoff surfaces perpendicular 
to each figure. It is apparent that the surfaces do not 
intersect any region where significant fluctuations are 
found. These appear to be the optimal locations and 
are used for the calculations. Similar reasoning is used 
in the minor axis plane. In a study of the sensitivity 
of the noise predictions to the location of the Kirchhoff 
surface, Lyrintzis 19 found that there is little difference 
between the results if the surface is placed at least one 
diameter away from the jet axis. 



Figure 8: The surface used to for the far field (Kirch- 
hoff) calculations. 



Figure 9: RMS of the mass flux perturbations in the 
major axis plane for the jet excited by a varicose mode. 
The values are normalized by the jet exit mass flux. 
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Figure 10: Far field SPL levels for the varicose excita- 
tion. 



Figure 11: Far field SPL levels for the flapping excita- 
tion. 

The observer locations are selected to match the mi- 
crophone stations used by Kinzie. 16 They consist of 
locations on an arc 25 D eq away from the center of the 
nozzle exit. The angular range is 15 to 50 degrees mea- 
sured from the jet axis. Overall sound pressure levels in 
the far field are presented in figures 10 and 11 for the 
varicose and flapping excitation cases and compared to 
the experimental data. The figures are a function of 
the angle the observer makes with the jet axis and are 
commonly referred to as directivity plots. 

For the locations close to the jet axis (i.e., small angle 
/?), the simulation overpredicts the experimental data 
for both the varicose and flapping cases. This is most 
likely due the higher amplitude instability waves (coher- 
ent structures) found in the simulated jet. Recall that 
the simulation predicts a peak instability wave ampli- 
tude that is approximately twice that of the experimen- 
tal jet. If the peak amplitude in the experimental jet 
is scaled by a factor of two, and one assumes that the 
noise is dominated by turbulent mixing noise, then the 
radiated pressure amplitude should be scaled by a factor 
of four (as deduced from the acoustic analogy 38 ). This 
results in a twelve decibel increase in the radiated noise. 
If the noise is dominated by Mach wave emission, then 
the instability wave analysis 39 predicts that the radi- 
ated pressure be scaled by a factor of two, which results 
in a 6 decibel increase in the far field noise. The peak 
differences between the simulation and experiment for 


the varicose and flapping excited jets are 7 and 12 dB, 
respectively. These differences suggest that the flapping 
jet has weak Mach wave radiation while the varicose jet 
noise is dominated by Mach wave radiation. 

The good correlation between the discrepancy in the 
noise source prediction and the discrepancy in the far 
field noise prediction leads one to conclude that the al- 
gorithm is capable of predicting the far field noise ra- 
diation for a given source amplitude. This argument 
is strengthened by the good agreement in the trends of 
the simulated and the experimental jets. Considering 
first the comparison for the jet excited by a varicose 
mode (figure 10), the angle of peak noise radiation is 
correctly predicted to be near 25 degrees. The noise 
is more directional in the simulated jet, however, com- 
pared to the experimental one. Recall from the intro- 
duction that turbulent mixing noise is fairly directional. 
The shock cell structure in the experimental jet may 
explain why its noise is less directional. This structure 
leads to the additional shock associated noise genera- 
tion mechanisms. Since this type of noise has upstream 
propagating components, it will tend to make the noise 
less directional when the contributions of all of the noise 
sources are combined. Although there is also a shock 
cell structure in the simulated jet excited by a varicose 
mode, it is comparatively weak. 

The spectra of the far field noise in the major axis 
plane is shown in figure 12 for the simulated jet ex- 
cited by a varicose mode. Like the hot wire spectra 
presented earlier, screech tones are apparent in the ex- 
perimental data. Ignoring these tones, however, like the 
experimental jet, the simulated one does predict a fairly 
broad spectral peak centered near a Strouhal number of 
0 . 2 . 

5 Conclusions 

In this paper a general algorithm for the prediction of 
supersonic jet noise is presented. The algorithm con- 
sists of the numerical simulation of the noise sources 
and sound radiation to the acoustic near field. The 
time dependent near field solution is then passed on to 
a Kirchhoff formulation to determine the far field noise. 
The algorithm has been applied to a perfectly expanded, 
cold, supersonic rectangular jet problem. This geome- 
try and set of flow conditions provide two advantages. 
The first is that turbulent mixing should be the sole 
noise generation mechanism present. This type of noise 
is found in the lower frequencies of the spectrum and 
therefore reduced the grid resolution requirements. The 
second advantage is that these conditions correspond to 
those of an experiment that has been conducted recently 
at Penn State. Thus, a comparison with the experi- 
mental data has been possible. This comparison has 
revealed that the algorithm is capable of reproducing 
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Figure 12: Far field noise spectra in the major axis plane 
for the flapping excitation. The numbers to the right 
of each graph denotes angular location in degrees. The 
heavy lines are the experimental spectra. 

many of the flow and noise features found in the ex- 
perimental jet. One discrepancy is that the simulation 
overpredicts the amplitude of the instability waves by a 
factor of two. Thus, future work on this project should 
investigate the reasons for this discrepancy and devise 
methods to eliminate it. 
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